#setwd("/home/jar/zxjar/paper_helper/userdata/0/34")
library(readr)

library(gtsummary)
library(survey)
library(haven)
library(tableone)
library(plyr)
library(tidyr)
# 链接：https://dplyr.tidyverse.org/reference/mutate.html
library(dplyr)
#library(tidyverse)
library(arsenal)

setwd('F:/Rproject/test-first-r')

paper.data<- read_csv("fromdata.csv")

# 防止四分位数出现错误直接把要分析的四分位数把na的直接默认给0

paper.data[which(is.na(paper.data$DietaryLZ)),'DietaryLZ'] <-0
paper.data[which(is.na(paper.data$SupplementLZ)),'SupplementLZ'] <-0
paper.data[which(is.na(paper.data$TotalLZ)),'TotalLZ'] <-0


# 根据前端传进来的字段进行执行表一数据分析的维度
#analyze.variable <- c("Waist","AnimalFluency","CERADdelayrecall","CERAD1","CERAD2","CERAD3","DSST","PIR","Age","SDMVPSU","SDMVSTRA","WTDRD1","Sex","Agegroup","Race","Alqgroup","Totalcalories","CERADtotal","Smokegroup","educationattainment","DietaryLZ","SupplementLZ","TotalLZ")
# 去除掉没用的列 整合数据
#paper.data <- paper.data[, analyze.variable]

NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")

# table1 右侧

 
 #tab1 开始
 
tbl_svysummary(NHANES_design,  missing = 'no',
               include = c(Agegroup, Sex,  Race, BMIgroup, Alqgroup, Smokegroup, educationattainment,
                           Age, PIR, BMXBMI, Waist, Totalcalories, DietaryLZ, SupplementLZ, TotalLZ,
                           CERAD1, CERAD2, CERAD3, CERADtotal, CERADdelayrecall, AnimalFluency, DSST),
               statistic = list(all_continuous()  ~ "{mean}(SE={mean.std.error})  {median}({p25}, {p75})" # 或者常用的："{mean} ({sd})"
                                 # all_categorical() ~ "{n_weighted} ({p}%)"
                                )
               )%>%as_flex_table()%>%flextable::save_as_html(path = '2.8table1.html')


# 进行测试多维度测试
varindex <- 1
# 这里写不成自动替换 在java 进行拼接 直接替换

 
 #tab2 开始
 
DietaryLZ.quantile.res <- svyquantile(~ DietaryLZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
adasd <-DietaryLZ.quantile.res$DietaryLZ[,'quantile'] 
if(adasd[2]==adasd[3]||adasd[2]==adasd[4]||adasd[2]==adasd[5]||adasd[3]==adasd[4]||adasd[3]==adasd[5]||adasd[4]==adasd[5]){
  print("存在相等进行各个减掉对应值")
  adasd[2] <-  adasd[2]-0.003
  adasd[3] <-  adasd[3]-0.002
  adasd[4] <-  adasd[4]-0.001
} 
paper.data$DietaryLZ.quantile.var <- cut(paper.data$DietaryLZ,
                                          breaks = adasd,
                                          labels = c('Q1','Q2','Q3','Q4'))

paper.data$DietaryLZ.quantile.var[which(is.na(paper.data$DietaryLZ.quantile.var))] <- 'Q1'
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")
tbl_svysummary(NHANES_design, by = DietaryLZ.quantile.var, missing = 'no',
               include = c(DietaryLZ,AnimalFluency,CERADdelayrecall,CERAD1,CERAD2,CERAD3,CERADtotal,DSST))%>% 
               # label = list(DietaryLZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%  modify_spanning_header(all_stat_cols() ~ "**DietaryLZ**")%>%              
  add_p()%>%
  as_flex_table() %>% # 导出 html
  flextable::save_as_html(path = 'Table2_Result.html') 
SupplementLZ.quantile.res <- svyquantile(~ SupplementLZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
adasd <-SupplementLZ.quantile.res$SupplementLZ[,'quantile'] 
if(adasd[2]==adasd[3]||adasd[2]==adasd[4]||adasd[2]==adasd[5]||adasd[3]==adasd[4]||adasd[3]==adasd[5]||adasd[4]==adasd[5]){
  print("存在相等进行各个减掉对应值")
  adasd[2] <-  adasd[2]-0.003
  adasd[3] <-  adasd[3]-0.002
  adasd[4] <-  adasd[4]-0.001
} 
paper.data$SupplementLZ.quantile.var <- cut(paper.data$SupplementLZ,
                                          breaks = adasd,
                                          labels = c('Q1','Q2','Q3','Q4'))

paper.data$SupplementLZ.quantile.var[which(is.na(paper.data$SupplementLZ.quantile.var))] <- 'Q1'
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")
tbl_svysummary(NHANES_design, by = SupplementLZ.quantile.var, missing = 'no',
               include = c(DietaryLZ,AnimalFluency,CERADdelayrecall,CERAD1,CERAD2,CERAD3,CERADtotal,DSST))%>% 
               # label = list(SupplementLZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%  modify_spanning_header(all_stat_cols() ~ "**SupplementLZ**")%>%              
  add_p()%>%
  as_flex_table() %>% # 导出 html
  flextable::save_as_html(path = 'Table3_Result.html') 
TotalLZ.quantile.res <- svyquantile(~ TotalLZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
adasd <-TotalLZ.quantile.res$TotalLZ[,'quantile'] 
if(adasd[2]==adasd[3]||adasd[2]==adasd[4]||adasd[2]==adasd[5]||adasd[3]==adasd[4]||adasd[3]==adasd[5]||adasd[4]==adasd[5]){
  print("存在相等进行各个减掉对应值")
  adasd[2] <-  adasd[2]-0.003
  adasd[3] <-  adasd[3]-0.002
  adasd[4] <-  adasd[4]-0.001
} 
paper.data$TotalLZ.quantile.var <- cut(paper.data$TotalLZ,
                                          breaks = adasd,
                                          labels = c('Q1','Q2','Q3','Q4'))

paper.data$TotalLZ.quantile.var[which(is.na(paper.data$TotalLZ.quantile.var))] <- 'Q1'
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")
tbl_svysummary(NHANES_design, by = TotalLZ.quantile.var, missing = 'no',
               include = c(DietaryLZ,AnimalFluency,CERADdelayrecall,CERAD1,CERAD2,CERAD3,CERADtotal,DSST))%>% 
               # label = list(TotalLZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%  modify_spanning_header(all_stat_cols() ~ "**TotalLZ**")%>%              
  add_p()%>%
  as_flex_table() %>% # 导出 html
  flextable::save_as_html(path = 'Table4_Result.html') 




 
 #tab3 开始
 
TotalLZ.quantile.res <- svyquantile(~ TotalLZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
paper.data$TotalLZ.quantile.var <- cut(paper.data$TotalLZ,
                                       breaks = TotalLZ.quantile.res$TotalLZ[,'quantile'],
                                       labels = c('Q1', 'Q2', 'Q3', 'Q4'))

NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA,
                           nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")
# Dietary.LZ
dr.age.adjusted.tbl <- svyglm(CERADdelayrecall ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'))



dr.fully.adjusted.tbl <- svyglm(CERADdelayrecall ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'))

# SupplementLZ
supp.age.adjusted.tbl <- svyglm(CERADdelayrecall ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'))


supp.fully.adjusted.tbl <- svyglm(CERADdelayrecall ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'))

# Total L.Z gcr
Total.age.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'))


Total.fully.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr

Totallz.age.Quartile.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'))


Totallz.fully.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                     design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'))


# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERADdelayrecall<-tbl_stack(
  tbls = list(dr.age.adjusted.tbl, dr.fully.adjusted.tbl,
              supp.age.adjusted.tbl, supp.fully.adjusted.tbl,Total.age.adjusted.tbl,Total.fully.adjusted.tbl,Totallz.age.Quartile.adjusted.tbl,Totallz.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))



# CERAD1

dr.age.adjusted.tbl <- svyglm(CERAD1 ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'))


dr.fully.adjusted.tbl <- svyglm(CERAD1 ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'))

# SupplementLZ
supp.age.adjusted.tbl <- svyglm(CERAD1 ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'))


supp.fully.adjusted.tbl <- svyglm(CERAD1 ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                  design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'))

# Total L.Z gcr

Total.age.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'))


Total.fully.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr

Totallz1.age.Quartile.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'))



Totallz1.fully.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'))

# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD1col<-tbl_stack(
  tbls = list(dr.age.adjusted.tbl, dr.fully.adjusted.tbl,
              supp.age.adjusted.tbl, supp.fully.adjusted.tbl,Total.age.adjusted.tbl,Total.fully.adjusted.tbl,Totallz1.age.Quartile.adjusted.tbl,Totallz1.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))



# cerad2
dr2.age.adjusted.tbl <- svyglm(CERAD2 ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'))

dr2.fully.adjusted.tbl <- svyglm(CERAD2 ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                 design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'))
# SupplementLZ
supp2.age.adjusted.tbl <- svyglm(CERAD2 ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'))

supp2.fully.adjusted.tbl <- svyglm(CERAD2 ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'))
# Total L.Z gcr
Total2.age.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'))

Total2.fully.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                    design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'))

# Total L.Z Quartile gcr
Totallz2.age.Quartile.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'))

Totallz2.fully.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD2col<-tbl_stack(
  tbls = list(dr2.age.adjusted.tbl, dr2.fully.adjusted.tbl,
              supp2.age.adjusted.tbl, supp2.fully.adjusted.tbl,Total2.age.adjusted.tbl,Total2.fully.adjusted.tbl,Totallz2.age.Quartile.adjusted.tbl,Totallz2.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))

# CERAD 3

dr3.age.adjusted.tbl <- svyglm(CERAD3 ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'))
dr3.fully.adjusted.tbl <- svyglm(CERAD3 ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                 design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'))
# SupplementLZ
supp3.age.adjusted.tbl <- svyglm(CERAD3 ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'))

supp3.fully.adjusted.tbl <- svyglm(CERAD3 ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                   design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'))
# Total L.Z gcr
Total3.age.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'))

Total3.fully.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                    design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'))
# Total L.Z Quartile gcr
Totallz3.age.Quartile.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'))
Totallz3.fully.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'))

# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
CERAD3col<-tbl_stack(
  tbls = list(dr3.age.adjusted.tbl, dr3.fully.adjusted.tbl,
              supp3.age.adjusted.tbl, supp3.fully.adjusted.tbl,Total3.age.adjusted.tbl,Total3.fully.adjusted.tbl,Totallz3.age.Quartile.adjusted.tbl,Totallz3.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))



alllei<-tbl_merge(tbls = list(CERADdelayrecall,CERAD1col,CERAD2col,CERAD3col),tab_spanner = c('CERAD: Score DelayedRecall','CERAD: Trial 1 Score','CERAD: Trial 2 Score','CERAD: Trial 3 Score'))
alllei%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table3.html")


 
 #tab4 开始
 
# Dietary.LZ
# DietaryLZ

tb4dr.age.adjusted.tbl <- svyglm(CERADtotal ~ DietaryLZ + Age, design = NHANES_design)%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb4dr.fully.adjusted.tbl <- svyglm(CERADtotal ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                   design = NHANES_design) %>%
  tbl_regression( include = c(DietaryLZ),
                  label = list(DietaryLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# SupplementLZ
tb4supp.age.adjusted.tbl <- svyglm(CERADtotal ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4supp.fully.adjusted.tbl <- svyglm(CERADtotal ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                     design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
tb4Total.age.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
tb4Total.fully.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb4Totallz.age.Quartile.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )


tb4Totallz.fully.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                        design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
tb4CERAD.total<-tbl_stack(
  tbls = list(tb4dr.age.adjusted.tbl, tb4dr.fully.adjusted.tbl,
              tb4supp.age.adjusted.tbl, tb4supp.fully.adjusted.tbl,tb4Total.age.adjusted.tbl,tb4Total.fully.adjusted.tbl,tb4Totallz.age.Quartile.adjusted.tbl,tb4Totallz.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))

# Animal Fluency score
Animaldr.age.adjusted.tbl <- svyglm(AnimalFluency ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
Animaldr.fully.adjusted.tbl <- svyglm(AnimalFluency ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# SupplementLZ
Animalsupp.age.adjusted.tbl <- svyglm(AnimalFluency ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'))
Animalsupp.fully.adjusted.tbl <- svyglm(AnimalFluency ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup +
                                          Smokegroup + PIR + educationattainment,
                                        design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z gcr
AnimalTotal.age.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotal.fully.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup +
                                           Smokegroup + PIR + educationattainment,
                                         design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
AnimalTotallz1.age.Quartile.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
AnimalTotallz1.fully.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                            design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
Animal.Fluencytab<-tbl_stack(
  tbls = list(Animaldr.age.adjusted.tbl, Animaldr.fully.adjusted.tbl,
              Animalsupp.age.adjusted.tbl, Animalsupp.fully.adjusted.tbl,AnimalTotal.age.adjusted.tbl,AnimalTotal.fully.adjusted.tbl,AnimalTotallz1.age.Quartile.adjusted.tbl,AnimalTotallz1.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))


# Digit Symbol Score
#Dietary L and Z
dss.age.adjusted.tbl <- svyglm(DSST ~ DietaryLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dss.fully.adjusted.tbl <- svyglm(DSST ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                 design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# SupplementLZ
dsssupp2.age.adjusted.tbl <- svyglm(DSST ~ SupplementLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dsssupp2.fully.adjusted.tbl <- svyglm(DSST ~ SupplementLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                      design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(SupplementLZ),
                 label = list(SupplementLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z gcr
dssTotal2.age.adjusted.tbl <- svyglm(DSST ~ TotalLZ + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
dssTotal2.fully.adjusted.tbl <- svyglm(DSST ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                       design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

# Total L.Z Quartile gcr
dssTotallz2.age.Quartile.adjusted.tbl <- svyglm(DSST ~ TotalLZ.quantile.var + Age, design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Age-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )

dssTotallz2.fully.adjusted.tbl <- svyglm(DSST ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,
                                         design = NHANES_design) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Fully-adjusted'),pvalue_fun = function(x) style_pvalue(x, digits = 2))%>%
  modify_header(
    std.error = "**SE**"
  )
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）
DSSTcol<-tbl_stack(
  tbls = list(dss.age.adjusted.tbl, dss.fully.adjusted.tbl,
              dsssupp2.age.adjusted.tbl, dsssupp2.fully.adjusted.tbl,dssTotal2.age.adjusted.tbl,dssTotal2.fully.adjusted.tbl,dssTotallz2.age.Quartile.adjusted.tbl,dssTotallz2.fully.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ", "SupplementLZ", "SupplementLZ","TotalLZ", "TotalLZ","Quartile of total L and Z","Quartile of total L and Z"))
# regression 的结果一般采用纵向拼接（stack）而不是横向拼接（merger）

alllei4<-tbl_merge(tbls = list(tb4CERAD.total,Animal.Fluencytab,DSSTcol),tab_spanner = c('CERAD: Total score','Animal Fluency score','Digit Symbol Score'))
alllei4%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table4.html")


 
 #tab5 开始
 
options( survey.lonely.psu = "adjust" )

# Dietary L and Z
CERAtotal.White.adjusted.tbl <- svyglm(CERADdelayrecall ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERADdelayrecall ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERADdelayrecall ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5sdr <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

#CERAD: Trial 1 Score 
CERAtotal.White.adjusted.tbl <- svyglm(CERAD1 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD1 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD1 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad1 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 


#CERAD: Trial 2 Score
CERAtotal.White.adjusted.tbl <- svyglm(CERAD2 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD2 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD2 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad2 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

#CERAD: Trial 3 Score
CERAtotal.White.adjusted.tbl <- svyglm(CERAD3 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

CERAtotal.black.adjusted.tbl <- svyglm(CERAD3 ~ DietaryLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Total L and Z
CERAtotalLZ.White.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalLZ.black.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
#Quartile of total L and Z
CERAtotalQUAN.White.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment,design = subset(NHANES_design, Race == 'Non-Hispanic White')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
CERAtotalQUAN.black.adjusted.tbl <- svyglm(CERAD3 ~ TotalLZ.quantile.var+ Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment ,design = subset(NHANES_design, Race == 'Non-Hispanic Black')) %>%
  tbl_regression(exponentiate = TRUE, include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
table5cerad3 <- tbl_stack(
  tbls = list(CERAtotal.White.adjusted.tbl, CERAtotal.black.adjusted.tbl,CERAtotalLZ.White.adjusted.tbl,CERAtotalLZ.black.adjusted.tbl,CERAtotalQUAN.White.adjusted.tbl,CERAtotalQUAN.black.adjusted.tbl),
  group_header = c("DietaryLZ", "DietaryLZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z")) 

alllei5<-tbl_merge(tbls = list(table5sdr,table5cerad1,table5cerad2,table5cerad3),tab_spanner = c('CERAD: Score Delayed Recall','CERAD: Trial 1 Score','CERAD: Trial 2 Score','CERAD: Trial 3 Score'))
alllei5%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table5.html")


 
 #tab6 开始
 
tb6drwhite.age.adjusted.tbl <- svyglm(CERADtotal ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6drblack.age.adjusted.tbl <- svyglm(CERADtotal ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L and Z
tb6talzwhite.age.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6talzblack.age.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb6Quartilewhite.age.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6Quartileblack.age.adjusted.tbl <- svyglm(CERADtotal ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
table6cts <- tbl_stack(
  tbls = list(tb6drwhite.age.adjusted.tbl, tb6drblack.age.adjusted.tbl,tb6talzwhite.age.adjusted.tbl,tb6talzblack.age.adjusted.tbl,tb6Quartilewhite.age.adjusted.tbl,tb6Quartileblack.age.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z"))
# Animal Fluency score

tb6drwhite.age.adjusted.tbl <- svyglm(AnimalFluency ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6drblack.age.adjusted.tbl <- svyglm(AnimalFluency ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L and Z
tb6talzwhite.age.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6talzblack.age.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb6Quartilewhite.age.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6Quartileblack.age.adjusted.tbl <- svyglm(AnimalFluency ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
table6af <- tbl_stack(
  tbls = list(tb6drwhite.age.adjusted.tbl, tb6drblack.age.adjusted.tbl,tb6talzwhite.age.adjusted.tbl,tb6talzblack.age.adjusted.tbl,tb6Quartilewhite.age.adjusted.tbl,tb6Quartileblack.age.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z"))
# Digit Symbol Score
tb6drwhite.age.adjusted.tbl <- svyglm(DSST ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6drblack.age.adjusted.tbl <- svyglm(DSST ~ DietaryLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(DietaryLZ),
                 label = list(DietaryLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L and Z
tb6talzwhite.age.adjusted.tbl <- svyglm(DSST ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6talzblack.age.adjusted.tbl <- svyglm(DSST ~ TotalLZ + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ),
                 label = list(TotalLZ ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
# Total L.Z Quartile gcr
tb6Quartilewhite.age.adjusted.tbl <- svyglm(DSST ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic White'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'White'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
tb6Quartileblack.age.adjusted.tbl <- svyglm(DSST ~ TotalLZ.quantile.var + Age + Sex + BMXBMI + Alqgroup + Smokegroup + PIR + educationattainment, design = subset(NHANES_design, Race == 'Non-Hispanic Black'))%>%
  tbl_regression(include = c(TotalLZ.quantile.var),
                 label = list(TotalLZ.quantile.var ~ 'Black'),
                 pvalue_fun = function(x) style_pvalue(x, digits = 2)
  )%>%
  modify_header(
    std.error = "**SE**"
  )
table6DSST <- tbl_stack(
  tbls = list(tb6drwhite.age.adjusted.tbl, tb6drblack.age.adjusted.tbl,tb6talzwhite.age.adjusted.tbl,tb6talzblack.age.adjusted.tbl,tb6Quartilewhite.age.adjusted.tbl,tb6Quartileblack.age.adjusted.tbl),
  group_header = c("Dietary.LZ", "Dietary.LZ","Total L and Z","Total L and Z","Quartile of total L and Z","Quartile of total L and Z"))


alllei6<-tbl_merge(tbls = list(table6cts,table6af,table6DSST),tab_spanner = c('CERAD: Total score','Animal Fluency score','Digit Symbol Score'))
alllei6%>%as_flex_table()%>%flextable::save_as_html(path = "2.8table6.html")




